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Computations of Flow over a Hump Model Using Higher Order 
Method with Turbulence Modeling 


P.Balakumar' 

NASA Langley Research Center, Hampton, VA 23681 

Turbulent separated flow over a two-dimensional hump is computed by solving the 
RANS equations with k-to (SST) turbulence model for the baseline, steady suction 
and oscillatory blowing/suction flow control cases. The flow equations and the 
turbulent model equations are solved using a fifth-order accurate weighted 
essentially nonoscillatory (WENO) scheme for space discretization and a third 
order, total variation diminishing (TVD) Runge-Kutta scheme for time integration. 
Qualitatively the computed pressure distributions exhibit the same behavior as 
those observed in the experiments. The computed separation regions are much 
longer than those observed experimentally. However, the percentage reduction in 
the separation region in the steady suction case is closer to what was measured in 
the experiment. The computations did not predict the expected reduction in the 
separation length in the oscillatory case. The predicted turbulent quantities are two 
to three times smaller than the measured values pointing towards the deficiencies in 
the existing turbulent models when they are applied to strong steady /unsteady 
separated flows. 


I. Introduction 

Computing unsteady separated turbulent flow’s accurately and efficiently is currently a 
challenging problem in CFD research. Most of the computational codes at present use implicit second 
order methods with one or two equation turbulence models with multi grid convergence acceleration 
techniques. Higher order methods are successfully used in direct numerical simulation (DNS) of stability 
and transition of shear flows, of turbulent shear flows at low Reynolds numbers and in Large Eddy 
Simulations (LES) of shear flows. Time integrations are performed using fourth order Runge-Kutta type 
algorithms and space discretizations are performed using fourth or higher order central, upwind, compact 
or spectral methods with selected higher order filtering when necessary. High operating Reynolds 
numbers and complex configurations hinder applying these schemes to engineering calculations and in 
practice modeled Reynolds Averaged Navier-Stokes (RANS) equations are solved to obtain engineering 
quantities. The current status and the issues related to higher order methods are discussed in Ref. 1 . 

In the finite difference formulation, higher order spatial discretization requires larger stencil size 
and smooth grids. This makes simple extension of the higher order methods to complex geometries 
difficult. Block structured or overset grid systems w’ith interpolation techniques near the interfaces are 
required to overcome these difficulties. Among the higher-order schemes, compact schemes 1 2 and upwind- 
based schemes’ are pursued as viable algorithms to develop higher order codes. Compact schemes are 
global in character and in general have symmetric coefficients, hence have low dispersion and dissipation 
properties and provide high resolution with minimum number of grid points. However, due to the low 
dissipation, high frequency oscillations are generated and artificial damping or filtering are required to 
remove these instabilities. Another issue with these schemes is that across flow discontinuities and near 
under resolved large gradients they exhibit the classical Gibbs phenomena and the schemes have to be 
modified to upwind biased flux splitting type schemes and this may be difficult to implement in a general 
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setting. Compact schemes have been applied to a range of problems 4, 5 ' 6 . Higher order upwind based 
methods are essentially non-oscillatory (ENO) type methods that are designed to suppress the oscillations 
near discontinuities and larger gradient regions. They are robust and can be applied to a wide range of 
Mach number flows without any modifications. However due to the upwind biased stencil and flux 
splitting they are more dispersive and dissipative compared to compact schemes. 

The objectives of this work are to apply a higher order accurate finite difference code to complex 
turbulent flows and to validate the computations against the experimental results. Three dimensional 
compressible unsteady Navier-Stokes equations are solved using a fifth order spatially accurate weighted 
essentially non-oscillatory (WENO) scheme and third-order accurate total variation diminishing (TVD) 
Runge-Kutta integration scheme. This is an upwind biased conservative scheme designed to suppress the 
oscillations near the discontinuities and large gradients. This method has been successfully applied to 
high Reynolds number laminar flows to investigate shock boundary layer interactions 7 , the stability and 
transition of supersonic boundary layers induced by roughness and acoustic waves 8 , and the screech tone 
generated in an imperfectly expanded supersonic jet 9 . Further development of this code is to extend this 
method to include turbulence modeling, and advanced modeling like DES and LES capabilities where 
low numerical dissipation is advantageous in discerning the effects of the sub-grid scale models from that 
introduced by numerical dissipation. In this paper, computations are performed for a turbulent flow over 
a two-dimensional hump model for the baseline, steady suction and oscillatory (suction/blowing) 
separation control cases by solving the RANS equations with k-to (SST) model using a higher-order finite 
difference code. A careful experiment has been performed for this configuration at NASA Langley 
Research center and the detailed experimental results are presented in the CFD validation workshop 
reports 10 and in Ref. 1 1. The model and the experiments are similar to the work performed by Seifert and 
Pack 12 . 


II. Experimental set up 

The experimental set up consists of a Glauret-Goldschmied type body mounted on a splitter plate. 
The chord length of the model is c = 16.53 inches and the height is 2.11 inches. The leading and the 
trailing edges are smoothly flared with the splitter plate and the height from the splitter plate to the upper 
wind tunnel wall is 15.03 inches. The splitter plate is 0.5 inches thick and extends 76.18 inches upstream 
and 44.45 inches downstream from the leading and trailing edges of the model. Tripping is applied at the 
leading edge of the splitter plate to obtain a fully developed turbulent boundary layer over the model. 
Steady suction and unsteady zero mass flux synthetic jet controls are applied through the slot opening 
located between 10.81-10.88 inches from the leading edge of the model. The model has two end-plates at 
both sides of the hump. The distance between the endplates is approximately 23 inches. The experiments 
are performed for a wide range of flow parameters such as Reynolds number, Mach number, different 
steady suction rates, different unsteady frequencies, and unsteady flow rates. An elaborate data 
acquisition system including streamwise and spanwise surface pressure measurements, pitot tube and hot 
wire measurements, 2-D and 3-D PIV and LDV and skin friction measurements were used to gather very 
detailed steady and unsteady flow field data. The experiments showed that the flow field is nominally 
two-dimensional in the center part of the model. The computations are performed for a two-dimensional 
model. 


III. Computational set up 

The computational domain extends from x/c = -10.0 to 4.0 in the streamwise direction and 
extends from the splitter plate to the upper tunnel wall in the normal direction. The leading edge of the 
splitter plate is modeled as a super ellipse with 0.25 inches half-thickness and an aspect ratio of 2. The 
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leading edge of the splitter plate is located at x/c = -5.9 and the leading edge of the hump model is located 
at x/c = 0.0. The length of the splitter plate is selected to match the measured velocity profiles at x/c - - 
2.1. A schematic diagram of the computational model is shown in Fig. 1. The free-stream Mach number 
is 0. 1 and the chord Reynolds number is 936.000. In the steady suction case, the suction flow rate is 23 
cfm through the suction slot across the span of 23 inches, and in the oscillatory' case the frequency of 
actuation is 137 Hz. and the maximum exit velocity is 27 m/'s. 


0.8 ; 


Inviscid wall 


0.6 | 

y/c ! inflow 


Outflow 




0.2 

0 


cn; Suction slot 

Elliptical x / c = o 65 

Symmetry leading edge Viscous walL~ 


K 


10 -8 -6 -4 x/ c -2 0 2 4 

Figure 1. Schematic diagram of the hump model. 

IV. Numerical Method 


The governing equations, the flow equations and the turbulent equations, are solved using the 5 th 
order accurate weighted essentially non-oscillatory (WENO) scheme for space discretization and using 
explicit third order total-variation-diminishing (TVD) Runge-Kutta scheme for time integration. The 
WENO and the TVD methods and the formulas are explained in Ref. 13 and the application of ENO 
method to N-S equations is given in Ref. 14. The solution method implemented in this computation is 
described in detail in Ref. 7. Standard k-co (SST) model is used and the equations and the model 
coefficients are described in Refs. 15-18. In the following presentation, (x, y) are the Cartesian 
coordinates, (u, v) are the velocity components, p is the density, and p is the pressure. E is the total energy 
given by 


r- li +Y~ 

E=e + — - — ,e = cJ.p = pRT. (1) 

Here e is the internal energy, T is the temperature, c y is the specific heat at constant volume and R is the 
gas constant. The viscosity (p ) is computed using Sutherland's law and the Prandtl number is taken as a 
constant value of 0.72. For turbulent quantities k is the turbulent kinetic energy, w is the dissipation rate 
and p T is the turbulent eddy viscosity. In the following subscript 'vc' denotes the values at the wall and the 
subscript denotes the free stream values. 

A C-Type grid is used around the splitter plate and a rectangular grid is added upstream of the 
leading edge as shown in Fig. 2. The rectangular grid overlaps the C grid and 5 th order central 
interpolation is used to transfer the flow variables from one grid to the other at the boundaries. A grid 
size of (651*151) is used around the splitter plate and the hump and a grid size of (101*51 ) is used in the 
rectangular region. 
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Figure 2: Overlapping grid near the leading edge of the splitter plate. 


A. Boundary conditions 

The following boundary conditions are implemented at different boundaries: 


1 . 


2 . 


3. 

4. 

5. 

6 . 


7 . 


At the upper wall inviscid conditions are applied. 
dp du dE 

-^ = — = — = v = 0. (2) 

dy dy dy 

At the lower wall viscous conditions are used. 
u = v = 0, 

'T' T 1 r T~' 13) 

1 w free stream ^ 

and p is computed from the continuity equations. 

From the leading edge of the splitter plate to the inflow boundary, symmetric conditions are used. 

At the inflow boundary stagnation pressure, one Riemann variable and normal velocity v=0 are 
prescribed and the second Riemann variable is solved to obtain the other flow quantity. 

At the outflow boundary the pressure is specified to obtain the required Mach number and 
characteristic-type boundary conditions are implemented similar to that described in Ref. 19 to 
obtain other flow variables. 

In the suction case, boundary conditions are applied on the surface of the hump across the suction 
slot. The suction slot extends from x/c = 0.654 to 0.658 and across the slot normal mass flow rate 
is specified to match the experimental value of 27.13 ft 3 /minute through the 23 inches span slot. 
A suction distribution of the form 


(P V )„ = /max 


(a end ^ start ) / 


(4) 


is used to get the required mass flow rate. Here f max is the maximum amplitude of the suction 
distribution selected to match the experimental total suction rate. Other forms have been tried 
and all of them yield the same results for a fixed total suction rate. The other flow quantities are 
obtained by extrapolation from inside the domain. 

In the oscillatory control case, boundary conditions similar to the steady suction case are 
implemented. A blowing/suction distribution in the form 


<T)„ = V max 


sirT 


' K{ X ~ X s,ar , ) ' 
\ (* en d ~ X start ) ) 


sin(o)/) 


(5) 


is prescribed for the velocity component normal to the slot. Here v max is the maximum normal 
oscillating velocity selected to match the experimental value. During the blowing cycle the 
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tangential velocity is obtained by assuming that the jet enters the flow domain at an angle of 30 
degrees to the surface of the hump, and during the suction cycle the other flow quantities are 
obtained by extrapolation from inside the domain. The maximum normal velocity and the 
frequency of oscillations are selected to match the experimental conditions of maximum exit 
velocity of 27 m/s and the oscillation frequency of 137 Hz. 


The following boundary conditions are implemented for the turbulent quantities at different boundaries. 


1 . 


2 . 

3. 

4. 


At the inflow boundary, small values are prescribed for k and Lij . 



= . 009 . 


( 6 ) 


At the outflow boundary', k and to are solved for from the governing equations. Higher order 
extrapolation condition is also tried and it gives the same results. 

At the lower viscous wall, k = 0 condition is used and the exact boundary condition is derived for 
to as described below. 

In the steady suction and oscillatory' control cases, across the suction slot, linear extrapolation is 
used to obtain the turbulent quantities on the surface. 


Since the variable to becomes singular near a viscous wall, in practice a large approximate value is 
prescribed at the wall. The turbulent model equations and the model coefficients are given in Ref. 1 8. 

60 


all ~ ' 


pM 


( 7 ) 


where d K is the distance to the first grid point from the wall. This is an approximate boundary condition 
and when it is implemented in the higher order scheme, oscillations and convergence problems are 

encountered and the following exact condition is derived for co wall . By realizing that a type 

singularity for the variable to arises because of the balance between the dissipation term and the viscous 
diffusion term in the to equation, the singularity is removed by rewriting the variable to as 

C 

to = —to v ( 8 ) 

y'n 

where y„ is the normal distance to the wall, C is a constant and to, is the new variable which is now 
regular near the wall. When this is substituted into the to equation the following equation is obtained for 
co h which is similar to the to equation except for the source term. 
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The value of co, at the wall becomes 


co, c?x / 


v / 


CO, 


1 wall 


6 K, dy n 

PP w \dx 


\ 


j) w 


( 9 ) 


(10) 


Here the variable co is nondimensionalised by 

* 

co 


co = 


K 2 ’ 


and C = 


Re 2 


(ID 


V v o 


Hence the procedure is to use the co, equation for the first few points near the wall and switch to the co 
equation away from the wall. In these computations, the co, equation is solved for the first ten points near 
the wall. Figure 3 shows the distribution of k, co and co, near the wall and it is seen that this technique 
resolves the viscous layer smoothly even though co is infinite at the wall. 



Figure 3: Variation of k, co and co I near the wall at xlc = 0.4. 
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V. Results 


Computations are performed for the baseline, the steady suction and the zero-mass oscillatory 
control cases. Figure 4 shows the comparison between the computed and the measured w-velocity profiles 
and the turbulent stress component (uu) profiles at x/c = -2.10. This is the furthest upstream location 
where the measurements have been made. The 'length of the splitter plate in the computational model is 
selected to match the computed w-velocity profiles with the experiment at this station. It is seen that the 
computed w-velocity profile agrees well with the measured w-velocity profile, but the turbulent quantity is 
under predicted by a factor of two in the computations. 
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Figure 4. Comparison of velocity and turbulent quantity profiles at x/c = -2.1 
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Figure 5: Contours of u velocity near the leading edge 
and in the entire computational domain. 


A. Baseline and steady suction control cases 

Figure 5 shows the contours of the w velocity near the leading edge region and for the entire 
computational domain. Near the leading edge region the flow separates and forms a small separation 
bubble. The smoothness of the solution in the overlap region shows the applicability of the interpolation 
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technique in the overlap region in this higher order formulation. Figure 6 shows the computed and the 
measured Cp distribution for the baseline and the steady suction cases. Qualitatively the same phenomena 
are observed in both cases and quantitatively there are some discrepancies in the pressure distributions in 
the separation region. 



x/c 

Figure 6. Cp distribution for the baseline and the steady suction control case. 

The flow decelerates as it approaches the leading edge of the hump and then steeply accelerates 
up to x/c = 0.20. In the baseline case, the acceleration is weakened after x/c = 0.20 and the pressure 
coefficient peaks around x/c = 0.40. Downstream the flow steeply decelerates up to x/c = 0.65 and 
separates. The computed peak Cp is about 10% smaller than the experimental value. This is due to the 
blockage of the end plates. In the experiments, the non-dimensionalization is done using the measured 
values at x/c = -2.10 which is upstream of the end plates. Computational studies performed 10 , one using an 
altered top wall contour to approximately account for the end-plates blockage and another by solving the 
full three-dimensional problem including the end-plates, showed better agreement of the peak pressure 
with the experimental data and supported the conjecture that there exists a small blockage effect caused 
by the end plates. It is also observed that some differences exist in the Cp distribution between the 
computation and the experiment in the separation region and the separation bubble is longer in the 
computation. In the steady suction control case, the flow continues to accelerate up to the suction slot x/c 
= 0.65, and then steeply decelerates and separates at the slot. As expected the separation region becomes 
smaller compared to the baseline case and as in the no flow control case the separation region is longer in 
the computation compared to the experiment. 

Figure 7 shows the streamline pattern for the computations and the experiments for the baseline 
and the steady suction cases. Table 2. gives the separation point and the reattachment points obtained 
from the computations and the experiments for the base line and the suction cases. In the baseline case, 
the length of the separation bubble is about 0.42 in the experiments while it is 0.61 in the computations. 
When the flow is controlled using suction, the separation bubble lengths are 0.24 and 0.32 for the 
experimental and computational cases, respectively. 
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Figure 7. Comparison of separation region for the baseline and the steady suction control case 

between the computation and the experiment. 


Figures 8 and 9 show the comparison between the computed and measured mean u-velocity 
profiles at x/c = 0.80, 0.90. 1.20 and 1.30 for the baseline and the suction cases. Near the separation point, 
x/c= 0.80 and 0.90. the agreement between the computation and the experiment is good and the 
discrepancy increases near the reattachment region. This conforms to the earlier observation that the 
separation region is larger in the computation than in the experiment. Figures 10 and 1 1 show the 
comparison for the turbulent stress component (mv). The agreement is not good for the turbulent quantities 
and the computations under predict the turbulent quantities by a factor of two to three. This points to the 
deficiencies in the turbulent model for the separated flows. 


no flow no flow 



Figure 8. Comparison of mean velocity profiles for the baseline case. 
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Figure 9. Comparison of mean velocity profiles for the suction case. 
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Figure 10. Comparison of turbulent stresses profiles for the baseline case. 
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Figure 11. Comparison of turbulent stresses profiles for the suction case. 
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B. Oscillatory control 

After the baseline and steady suction control cases were computed time accurate computations 
w'ere performed for the zero-mass flux blowing and suction control cases. Since this is an explicit code, 
the numerical instability limits the time step to a small value. It takes about 10 6 time steps per cycle. The 
computations are continued until a periodic state in time is reached. Figure 12 shows the pressure as a 
function of phase angle (or time) for the last two cycles at a fixed location on the wall inside the separated 
region at x/c = 0.8. It takes about 5-6 cycles before a periodic state is reached. 



Figure 12. Pressure fluctuations as a function phase angle (or time) 
at a fixed point on the wall xlc = 0.8. 

Figure 13 shows the computed and the measured Cp distributions for the baseline and the 
oscillatory cases. The pressure distribution for the oscillatory case is the time averaged pressure for one 
cycle. Qualitatively the same phenomena are observed in both cases and quantitatively there are some 
discrepancies in the pressure distributions in the separation region. In the experiment, the pressure 
decreases immediately downstream of the slot, then increases up to the reattachment point. In the 
computations, it remains flat near the slot and then increases. As in other cases, the separation region is 
narrower in the experiment compared to the computations. 
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x/c 

Figure 13. Cp distribution for the baseline and the oscillatory control case. 

Figure 14 shows the mean streamline patterns in the separation zone obtained from the 
computation and the experiment. Clearly the separation region is smaller, about 20%, in the experiment 
compared to the calculations. Table 2 summarizes the experimental and the computed separation and 
reattachment points for all three cases. The reattachment point for the unsteady control case is about x/c = 
1.21 in the simulation and about 0.98 in the experiment. In the experiment, the attachment point moved 
from x/c = 1.10 for the baseline case to x/c = 0.98 for the oscillatory control case and in the computations 
the attachment point moves from x/c = 1.276 in the baseline case to x/c = 1.21 in the control case. Even 
though the trend is in the right direction, the reduction in the separation is more in the experiment. 




Figure 14. Comparison of separation region for the oscillatory control case. 
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Table 2 Separation and reattachment locations 



Separation 
Point xlc 

Reattachment 
Point xlc 

Baseline 

Computation 

0.662 

1.276 

Experiment 

0.680 

1.100 

Suction 

Computation 

0.682 

1.004 

Experiment 

0.700 

0.940 

Oscillatory 

Computation 

ujn | 


Experiment 

0.675 

0.98 


Figures 15 compares the unsteady flow fields obtained from the computations and the 
experiments at different phases of the control cycle. The figure shows the contours of the w-velocity and 
the instantaneous streamline patterns at different phase angles 0 = 0. 90. 180 and 270 degrees. Below each 
of the contour plots, the unsteady pressure coefficient, Cp-Cp mea „, at the corresponding phase angles are 
depicted. Here Cp mean is the time averaged pressure coefficient shown in Fig. 13. The phase angles, 0 = 0 
and 180 degrees, correspond to the start of the blowing and suction phases of the control and 0 = 90 and 
270 degrees correspond to maximum blowing and suction phases. Qualitatively, the computed flow field 
structure and the unsteady pressure variations are similar to that are observed in the experiments. At the 
start of the cycle 0 = 0, both in the experiments and computations, a small vortex forms near the slot and a 
large vortex precedes this small vortex downstream. Following the sequence, it is seen that these vortices 
become larger and move downstream with increasing phase angles. At the end of the blowing phase or the 
start of the suction of the suction phase <p = 180. there exist two co-rotating vortices in the separation 
zone. During the suction phase, the vortex closer to the slot becomes larger and the second vortex 
becomes weaker. The significant difference between the computation and the experiment is the 
appearance of a large vortex in the computation near the reattachment region. This vortex persists as an 
elongated vortex during the blowing phase. This implies that the suction phase is not as effective in the 
computations as in the experiments in pulling the separation region up towards the slot. It may be that 
there is not sufficient turbulent mixing in the reattachment region to bring the high momentum fluid 
towards the wall. This may be the reason for over predicting the separation region in the computations. 

VI. Conclusions 

Turbulent separated flow over a two-dimensional hump is computed by solving the RANS 
equations with k-a>( SST) turbulence model using a fifth-order accurate WENO method for the baseline, 
steady suction and zero mass oscillatory' control cases. Overlapping grids are used to compute the leading 
edge region. The singularity in the w equation near the wall is removed by solving the w equation in the 
transformed plane and this also removes the ambiguity in applying the wall boundary condition for co that 
exists now in practice. 

Qualitatively the computed pressure distributions exhibit the same behavior as those were 
observed in the experiments. The computed separation regions are much longer in the computations than 
in the experiments. Computations verified that steady suction could be used to significantly reduce the 
separated region in a flow. However, the computations predict a small reduction in separated region 
compared to the experiment. The predicted turbulent quantities are two to three times smaller than that are 
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measured and it points towards the deficiencies in the existing turbulent models when they are applied to 
strong separated steady and unsteady flows. 

Time accurate computations of the unsteady zero mass oscillatory control case revealed similar 
flow field structures and unsteady pressure variations in the separated region. However, the computations 
did not show the expected large percentage reduction in the separated region compared to the experiment. 
In the computations an additional vortex appears near the reattachment region and in the experiment this 
was not observed. This may be due to less turbulent mixing in this region. 
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Figure 15. Comparison of the flow field and the oscillatory pressure in the separation region for the 

oscillatory control. 
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